// based on Smith & Duller (1959) phenomenological model
// modified from matlab version: Geant4_Muon_histogram_generation_v2.m

/*
%   For non-standard energy and angular distributions Geant4 requires the 
%   utilization of the General Particle Source module which via the 
%   G4GeneralParticleSource class allows specification of user defined 
%   angular and energy distributions. The user defined histograms are 
%   specified in macro files using the following commands:
%
%   1.	/gps/particle mu+ 	      (specify particle type)
%   2.	/gps/ang/type user 	      (user defined histogram)
%   3.	/gps/hist/type theta 	  (zenith angle histogram)
%   4.	/gps/hist/point Bt Wt 	  (angular histogram values)
%   5.	/gps/ene/type Arb 	      (user defined histogram)
%   6.	/gps/hist/type arb 	      (point-wise energy spectrum)
%   7.	/gps/hist/point Eh Hh 	  (energy spectrum values)
%   8.	/gps/hist/inter Lin 	  (interpolation scheme: Linear)
%   9.	/run/beamOn 1000	      (number of particles)
%
%   where a short explanation of each command appears in parentheses. 
%   The histograms represent differential functions and must be included 
%   one bin at a time. Angular histogram is described using the bin upper 
%   boundary and the area of the bin. Energy spectrum (point-wise) is 
%   described using the bin center and the height of the bin. The first 
%   value of each histogram must be the lower boundary of the bin and a 
%   dummy value that is not used.
%
%   Instructions: Run the file. On the command line set the minimum and 
%   maximum muon energy and the minimum and maximum zenith angle. 
%   Output:
%
%   1. The output "Muon_table" matrix contains the angles and energies 
%      of the sampled muons
%   2. The "User_defined_hist_energy" matrix contains the point-wise energy
%      (MeV) spectrum for the Geant4 macro file. Just copy paste to your Geant4 
%      macro file.
%   3. The "User_defined_hist_angle" matrix contains the zenith
%      angle (radians) histogram for the Geant4 macro file. Just copy paste to your 
%      Geant4 macro file.
%   4. The "User_defined_hist_angle_corrected" matrix contains the zenith
%      angle (radians) - corrected for the solid angle effect- histogram for 
%      the Geant4 macro file. Just copy paste to your Geant4 macro file.
%
%   This file is free for use. More details, examples and validation results 
%   can be found on the journal papers:
%   
%   S. Chatzidakis, S. Chrysikopoulou, L.H. Tsoukalas (2015)
%   "Developing a cosmic ray muon sampling capability for muon tomography
%   and monitoring applications", Nuclear Instruments and Methods in 
%   Phyics Research Section A, Vol. 804, pp. 33-42.
%
%   S. Chatzidakis, L.H. Tsoukalas (2015)
%   "A Geant4-MATLAB muon generator for Monte-Carlo simulations", 
%   URL: https://engineering.purdue.edu/~aisl/Stylianos_Publications.html
%
%   Users are kindly requested to cite the above journal papers in their
%   publications. Please comment. Thank you!
%
*/

void gen()
{
    double Emin = 1;   //minimum muon energy (GeV)
    double Emax = 100; //maximum muon energy (GeV)

    double Theta_min = 0.;  //minimum zenith angle(raneg 0~90) in degrees
    double Theta_max = 70.; //maximum zenith angle

    int N = 100000; ///number of muons to be sampled

    //constant
    double Alpha = 0.002382; // constant A
    double lambda = 120;     // absorption mean free path 120 g/cm2
    double kappa = 2.645;    // exponent (-)
    double bp = 0.771;
    double jp = 148.16;    // correction factor (-); factor (GeV)
    double alpha = 0.0025; // muon energy loss in GeV/g/cm2
    double rho = 0.76;     // fraction of pion energy that is transferred to muon
    double y0 = 1000;      // atmoshperic depth g/cm2
    double bm = 0.8;
    double Bm = 1.041231831; // correction factor (-);
    double pi = 3.1415926;

    // Inverse trasform and Accept-Reject method
    // 1. Inverse trasform to sample muon zenith angle
    //    select muon angle (in radians) - corrected for solid angle effect - using inverse transform
    double theta1 = Theta_min * pi / 180;
    double theta2 = Theta_max * pi / 180;
    for (int i = 1; i < N; i++)
    {
        double tm = theta1 + (theta2 - theta1)*rand();
        double tm1 = 2 * tm / (pi);
        double theta_corr(i) = acos((1 - tm1) ^ (1 / 3));
    }
}